Fin most abundant taxa across sites:

sapply(c(biobakery_function, alpha, beta, taxa, norm, heat, varia), source, chdir = TRUE)

Data loading:

input <- list(metpahlan_ps = "~/Documents/GitHub/oral-microbiome/data/processed/metaphlan_ps.RDS",
              motus_ps = "~/Documents/GitHub/oral-microbiome/data/processed/motus/physeq.RDS",
              pathway_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_pw_ps.RDS",
              kegg_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_kegg_ps.RDS",
              l4c_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_l4c_ps.RDS",
              go_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_go_ps.RDS",
              infogo_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_infogo_ps.RDS",
              metac_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_metac_ps.RDS",
              pfam_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_pfam_ps.RDS")
purrr::map(input, 
           readRDS) -> objects
sample_data(objects$metpahlan_ps$physeq)$Site_Health <- paste0(sample_data(objects$metpahlan_ps$physeq)$Oral_Site,
                                                               "_",
                                                               sample_data(objects$metpahlan_ps$physeq)$Health_status)

Fig1 Identification of predominant streptococcus species at all sites

This figure is good, but with the approach we decided focusing solely on streptococcus species We should change it a bit so that it shows top 10 streptococcus species across all three sites

objects$metpahlan_ps$physeq -> tmp
# subset_samples(Health_status == "Healthy control") %>%
# subset_samples(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_tongue_strep
tmp %>%
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq::subset_taxa(Species %in% ma_saliva_tongue_strep) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "Subject",#treatment
                          facet_by = c("Oral_Site", "Health_status", "Subject"),
                          ntax = 100,
                          tax_aggregate = "Species",
                          tax_add = NULL) -> p
## Warning: There are only 16 taxa, showing all
p + facet_grid(. ~ Oral_Site + Health_status, scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 25, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 25,50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.1))  + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) + 
  theme(axis.text.y = element_text(angle = 0,  size = 8))

Identification of predominant streptococcus species at all sites

Fig1A: predominant streptococcus species across all samples in health

objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") -> tmp
# subset_samples(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_tongue_strep
tmp %>%
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq::subset_taxa(Species %in% ma_saliva_tongue_strep) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "Subject",#treatment
                          facet_by = c("Oral_Site", "Health_status", "Subject"),
                          ntax = 100,
                          tax_aggregate = "Species",
                          tax_add = NULL) -> p
## Warning: There are only 16 taxa, showing all
p + facet_grid(. ~ Oral_Site + Health_status, scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 25, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 25,50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.1))  + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) + 
  theme(axis.text.y = element_text(angle = 0,  size = 8))

### Fig1B: predominant streptococcus species saliva vs. tongue samples in health

objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") %>% 
  subset_samples(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_tongue_strep
tmp %>%
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq::subset_taxa(Species %in% ma_saliva_tongue_strep) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "Subject",#treatment
                          facet_by = c("Oral_Site", "Health_status", "Subject"),
                          ntax = 100,
                          tax_aggregate = "Species",
                          tax_add = NULL) -> p
## Warning: There are only 13 taxa, showing all
p + facet_grid(. ~ Oral_Site + Health_status, scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 25, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 25,50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.1))  + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) + 
  theme(axis.text.y = element_text(angle = 0,  size = 8))

ma_saliva_tongue_strep
##  [1] "Streptococcus_australis"     "Streptococcus_cristatus"    
##  [3] "Streptococcus_gordonii"      "Streptococcus_infantis"     
##  [5] "Streptococcus_mitis"         "Streptococcus_mutans"       
##  [7] "Streptococcus_oralis"        "Streptococcus_parasanguinis"
##  [9] "Streptococcus_salivarius"    "Streptococcus_sanguinis"    
## [11] "Streptococcus_sp_A12"        "Streptococcus_sp_F0442"     
## [13] "Streptococcus_sp_HMSC071D03"
ma_saliva_tongue_strep_sel <- c("Streptococcus_parasanguinis",
                                "Streptococcus_salivarius",
                                "Streptococcus_infantis")

Fig1C: predominant streptococcus species saliva vs. plaque in health

objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") %>% 
  subset_samples(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_plaque_strep
tmp %>%
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq::subset_taxa(Species %in% ma_saliva_plaque_strep) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "Subject",#treatment
                          facet_by = c("Oral_Site", "Health_status", "Subject"),
                          ntax = 100,
                          tax_aggregate = "Species",
                          tax_add = NULL) -> p
## Warning: There are only 13 taxa, showing all
p + facet_grid(. ~ Oral_Site + Health_status, scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 25, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 25,50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.1))  + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) + 
  theme(axis.text.y = element_text(angle = 0,  size = 8))

ma_saliva_plaque_strep
##  [1] "Streptococcus_anginosus_group"   "Streptococcus_australis"        
##  [3] "Streptococcus_cristatus"         "Streptococcus_gordonii"         
##  [5] "Streptococcus_infantis"          "Streptococcus_mitis"            
##  [7] "Streptococcus_mutans"            "Streptococcus_oralis"           
##  [9] "Streptococcus_parasanguinis"     "Streptococcus_salivarius"       
## [11] "Streptococcus_sanguinis"         "Streptococcus_sobrinus"         
## [13] "Streptococcus_sp_oral_taxon_056"
ma_saliva_plaque_strep_sel <- c("Streptococcus_oralis",
                                "Streptococcus_sanguinis",
                                "Streptococcus_gordonii")

Fig2 Overall bacterial activity of selected species

Load pathway abundance table and have a look:

objects$pathway_ps$physeq  %>%
  humann2_species_contribution(meta_data_var = c("Subject", "sample_Sample" ,"Type", "Oral_Site", "Health_status")) %>%
  separate(Feature, into = c("code", "name"), sep = ": ", remove = FALSE) %>%
  dplyr::filter(!is.na(name))  -> pwy
## Warning in speedyseq::psmelt(.): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 720 rows [1029,
## 1903, 1961, 2338, 2461, 2854, 3562, 3891, 4578, 5649, 5690, 5849, 5854, 6607,
## 6843, 7028, 7077, 7309, 7819, 8576, ...].
pwy %>%
  distinct(name, .keep_all = TRUE) %>%
  DT::datatable()

Fig2A: box plot species saliva vs. tongue

ma_saliva_tongue_strep_sel
## [1] "Streptococcus_parasanguinis" "Streptococcus_salivarius"   
## [3] "Streptococcus_infantis"
pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = ". ~ .",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot + ggpubr::rotate()

plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Oral_Site,
                        group.by = c("Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

Fig2B: log DNA/Log RNA plots of species saliva vs. tongue

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_plot(x_plot = "log10(DNA)", 
                       y_plot = "log10(RNA)", 
                       color = "Species", 
                       fill = "Species",
                       shape = "Genus",
                       group = "Species",
                       filter_species = ma_saliva_tongue_strep_sel,
                       facet_formula = "Oral_Site ~ .") -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot 

plot$plot$data %>%
  select(-id, -sample_Sample, Health_status) %>%
  group_by(Subject, Oral_Site) %>%
  add_count(name = "count_per") %>% 
  slice_max(order_by = RNA, n = 3) %>%
DT::datatable()
plot$plot$data %>%
  select(-id, -sample_Sample, Health_status) %>%
  group_by(Subject, Oral_Site) %>%
  add_count(name = "count_per") %>% 
  slice_max(order_by = RNA_DNA, n = 3) %>%
DT::datatable()

Fig2C: box plot species saliva vs. plaque

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_plaque_strep_sel,
                             facet_formula = ". ~ .",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot + ggpubr::rotate()

plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Oral_Site,
                        group.by = c("Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

Fig2D: log DNA/Log RNA plots of species saliva vs. plaque

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  humann2_RNA_DNA_plot(x_plot = "log10(DNA)", 
                       y_plot = "log10(RNA)", 
                       color = "Species", 
                       fill = "Species",
                       shape = "Genus",
                       group = "Species",
                       filter_species = ma_saliva_plaque_strep_sel,
                       facet_formula = "Oral_Site ~ .") -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot 

plot$plot$data %>%
  select(-id, -sample_Sample, Health_status) %>%
  group_by(Subject, Oral_Site) %>%
  add_count(name = "count_per") %>% 
  slice_max(order_by = RNA, n = 3) %>%
DT::datatable()
plot$plot$data %>%
  select(-id, -sample_Sample, Health_status) %>%
  group_by(Subject, Oral_Site) %>%
  add_count(name = "count_per") %>% 
  slice_max(order_by = RNA_DNA, n = 3) %>%
DT::datatable()

S1A – pathway expression of selected streptococcus species in saliva vs. tongue in health

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

FigS1B: pathway expression of selected streptococcus species in saliva vs. plaque in health

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_plaque_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()

Fig 3. Pathway activity of selected species

Fig3A: significant pathway expression by candidate species saliva versus tongue in health

redo main plot + stats

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                       y_plot = "log2(RNA_DNA)",
                       color = "Oral_Site",
                       fill = "Oral_Site",
                       shape = NULL,
                       filter_species = "Streptococcus_infantis",
                       filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Health_status ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Oral_Site", 
                       fill = "Oral_Site",
                       shape = NULL,
                       filter_species = "Streptococcus_parasanguinis",
                       filter_feature = diff_signif %>% filter(Species == "Streptococcus_parasanguinis") %>% pull(Feature),
                       facet_formula = " . ~ Health_status ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Oral_Site", 
                       fill = "Oral_Site",
                       shape = NULL,
                       filter_species = "Streptococcus_salivarius",
                       filter_feature = diff_signif %>% filter(Species == "Streptococcus_salivarius") %>% pull(Feature),
                       facet_formula = " . ~ Health_status ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

Fig3B: significant pathway expression by candidate species saliva versus plaque in health

redo main plot + stats

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_plaque_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                       y_plot = "log2(RNA_DNA)",
                       color = "Oral_Site",
                       fill = "Oral_Site",
                       shape = NULL,
                       filter_species = "Streptococcus_oralis",
                       filter_feature = diff_signif %>% filter(Species == "Streptococcus_oralis") %>% pull(Feature),
                       facet_formula = " . ~ Health_status ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Oral_Site", 
                       fill = "Oral_Site",
                       shape = NULL,
                       filter_species = "Streptococcus_sanguinis",
                       filter_feature = diff_signif %>% filter(Species == "Streptococcus_sanguinis") %>% pull(Feature),
                       facet_formula = " . ~ Health_status ",
                       export_legend = TRUE) -> plot
  

plot$legend

plot$plot

# pwy %>%
#   dplyr::filter(Health_status == "Healthy control") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Oral_Site", 
#                        fill = "Oral_Site",
#                        shape = NULL,
#                        filter_species = "Streptococcus_gordonii",
#                        filter_feature = diff_signif %>% filter(Species == "Streptococcus_gordonii") %>% pull(Feature),
#                        facet_formula = " . ~ Health_status ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Fig 5. Overall bacterial activity of selected species in health versus diseases

Fig5A: box plot species in health and disease at each site

strep_sel <- c(ma_saliva_plaque_strep_sel, ma_saliva_tongue_strep_sel)
pwy %>%
  # dplyr::filter(Health_status == "Healthy control",
                # Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = strep_sel,
                             facet_formula = ". ~ Oral_Site",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot + ggpubr::rotate()

plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Health_status,
                        group.by = c("Oral_Site", "Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

Fig5B: log DNA/Log RNA plots of species in health vs. disease at each site

# pwy %>%
#   dplyr::filter(Health_status == "Healthy control",
#                 Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
#   humann2_RNA_DNA_plot(x_plot = "log10(DNA)", 
#                        y_plot = "log10(RNA)", 
#                        color = "Species", 
#                        fill = "Species",
#                        shape = "Genus",
#                        group = "Species",
#                        filter_species = strep_sel,
#                        facet_formula = "Oral_Site ~ .") -> plot

pwy %>%
   humann2_RNA_DNA_plot(x_plot = "log10(DNA)", 
                       y_plot = "log10(RNA)", 
                       color = "Species", 
                       fill = "Species",
                       shape = "Genus",
                       group = "Species",
                       filter_species = strep_sel,
                             facet_formula = "Health_status ~ Oral_Site") -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot$data %>%
  select(-id, -sample_Sample, Health_status) %>%
  group_by(Subject, Oral_Site) %>%
  add_count(name = "count_per") %>% 
  slice_max(order_by = RNA, n = 3) %>%
DT::datatable()
plot$plot$data %>%
  select(-id, -sample_Sample, Health_status) %>%
  group_by(Subject, Oral_Site) %>%
  add_count(name = "count_per") %>% 
  slice_max(order_by = RNA_DNA, n = 3) %>%
DT::datatable()

Fig 6 Pathway activity of selected species at each site in health versus disease

Fig6A-F: significant pathway expression by candidate species at each site in health and disease

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     # dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       # filter_species = "Streptococcus_infantis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species + Oral_Site ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

A - SALIVA

strep_sel <- c(ma_saliva_plaque_strep_sel, ma_saliva_tongue_strep_sel)
strep_sel
## [1] "Streptococcus_oralis"        "Streptococcus_sanguinis"    
## [3] "Streptococcus_gordonii"      "Streptococcus_parasanguinis"
## [5] "Streptococcus_salivarius"    "Streptococcus_infantis"

Streptococcus_infantis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_infantis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") -> diff

diff %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_infantis") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_parasanguinis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_parasanguinis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
   dplyr::filter(Species %in% "Streptococcus_parasanguinis") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       # filter_species = "Streptococcus_infantis",
                       filter_feature = diff_signif %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

Streptococcus_salivarius

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_salivarius",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
   dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       # filter_species = "Streptococcus_infantis",
                       filter_feature = diff_signif %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

Streptococcus_gordonii

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_gordonii",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_gordonii") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_sanguinis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_sanguinis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_sanguinis") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_oralis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_oralis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
   dplyr::filter(Species %in% "Streptococcus_oralis") %>%
     dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       # filter_species = "Streptococcus_infantis",
                       filter_feature = diff_signif %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

B - TONGUE_BIOFILM

strep_sel <- c(ma_saliva_tongue_strep_sel)
strep_sel
## [1] "Streptococcus_parasanguinis" "Streptococcus_salivarius"   
## [3] "Streptococcus_infantis"

Streptococcus_infantis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_infantis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_infantis") %>%
#      dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#   humann2_RNA_DNA_ratio_plot(x_plot = "name",
#                        y_plot = "log2(RNA_DNA)",
#                        color = "Health_status",
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif  %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
# 
# 
# plot$legend
# plot$plot

Streptococcus_parasanguinis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_parasanguinis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
   dplyr::filter(Species %in% "Streptococcus_parasanguinis") %>%
     dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                       y_plot = "log2(RNA_DNA)",
                       color = "Health_status",
                       fill = "Health_status",
                       shape = NULL,
                       # filter_species = "Streptococcus_infantis",
                       filter_feature = diff_signif %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

Streptococcus_salivarius

pwy %>%
   dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
     dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_salivarius",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
#      dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

C - SUBGINGIVAL_PLAQUE

strep_sel <- c(ma_saliva_plaque_strep_sel)
strep_sel
## [1] "Streptococcus_oralis"    "Streptococcus_sanguinis"
## [3] "Streptococcus_gordonii"

Streptococcus_oralis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_oralis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_infantis") %>%
#      dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#   humann2_RNA_DNA_ratio_plot(x_plot = "name",
#                        y_plot = "log2(RNA_DNA)",
#                        color = "Health_status",
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif  %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
# 
# 
# plot$legend
# plot$plot

Streptococcus_sanguinis

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_sanguinis",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_sanguinis") %>%
#      dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_gordonii

pwy %>%
   dplyr::filter(Species %in% strep_sel) %>%
     dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
   # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                       y_plot = "log2(RNA_DNA)", 
                       color = "Health_status", 
                       fill = "Health_status",
                       shape = NULL,
                       filter_species = "Streptococcus_gordonii",
                       # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                       facet_formula = " . ~ Species ",
                       export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
#      dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot
sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggConvexHull_0.1.0   ampvis2_2.6.4        reshape2_1.4.4      
##  [4] scales_1.1.1         phyloseq_1.34.0      forcats_0.5.1       
##  [7] stringr_1.4.0        dplyr_1.0.4          purrr_0.3.4         
## [10] readr_1.4.0          tidyr_1.1.2          tibble_3.1.0        
## [13] ggplot2_3.3.3        tidyverse_1.3.0.9000
## 
## loaded via a namespace (and not attached):
##   [1] ggbeeswarm_0.6.0     Rtsne_0.15           colorspace_2.0-0    
##   [4] ggsignif_0.6.1       rio_0.5.26           ellipsis_0.3.1      
##   [7] XVector_0.30.0       fs_1.5.0             rstudioapi_0.13     
##  [10] ggpubr_0.4.0         farver_2.1.0         fantaxtic_0.1.0     
##  [13] ggrepel_0.9.1        ggnet_0.1.0          DT_0.17             
##  [16] fansi_0.4.2          lubridate_1.7.10     xml2_1.3.2          
##  [19] codetools_0.2-18     splines_4.0.4        doParallel_1.0.16   
##  [22] knitr_1.31           ade4_1.7-16          jsonlite_1.7.2      
##  [25] broom_0.7.5          cluster_2.1.1        dbplyr_2.1.0        
##  [28] compiler_4.0.4       httr_1.4.2           backports_1.2.1     
##  [31] assertthat_0.2.1     Matrix_1.3-2         lazyeval_0.2.2      
##  [34] cli_2.3.1            htmltools_0.5.1.1    prettyunits_1.1.1   
##  [37] tools_4.0.4          igraph_1.2.6         gtable_0.3.0        
##  [40] glue_1.4.2           Rcpp_1.0.6           carData_3.0-4       
##  [43] Biobase_2.50.0       cellranger_1.1.0     jquerylib_0.1.3     
##  [46] vctrs_0.3.6          Biostrings_2.58.0    rhdf5filters_1.2.0  
##  [49] multtest_2.46.0      ape_5.4-1            nlme_3.1-152        
##  [52] iterators_1.0.13     crosstalk_1.1.1      xfun_0.21           
##  [55] network_1.16.1       openxlsx_4.2.3       rvest_0.3.6         
##  [58] lifecycle_1.0.0      rstatix_0.7.0        zlibbioc_1.36.0     
##  [61] MASS_7.3-53.1        hms_1.0.0            parallel_4.0.4      
##  [64] biomformat_1.7.0     rhdf5_2.34.0         RColorBrewer_1.1-2  
##  [67] curl_4.3             yaml_2.2.1           speedyseq_0.5.3.9001
##  [70] sass_0.3.1           stringi_1.5.3        highr_0.8           
##  [73] S4Vectors_0.28.1     foreach_1.5.1        permute_0.9-5       
##  [76] BiocGenerics_0.36.0  zip_2.1.1            rlang_0.4.10        
##  [79] pkgconfig_2.0.3      evaluate_0.14        lattice_0.20-41     
##  [82] Rhdf5lib_1.12.1      labeling_0.4.2       patchwork_1.1.1     
##  [85] htmlwidgets_1.5.3    cowplot_1.1.1        tidyselect_1.1.0    
##  [88] plyr_1.8.6           magrittr_2.0.1       R6_2.5.0            
##  [91] IRanges_2.24.1       generics_0.1.0       DBI_1.1.1           
##  [94] foreign_0.8-81       pillar_1.5.0         haven_2.3.1         
##  [97] withr_2.4.1          mgcv_1.8-34          abind_1.4-5         
## [100] survival_3.2-7       car_3.0-10           modelr_0.1.8        
## [103] crayon_1.4.1         utf8_1.1.4           microbiome_1.10.0   
## [106] plotly_4.9.3         rmarkdown_2.7        progress_1.2.2      
## [109] grid_4.0.4           readxl_1.3.1         data.table_1.14.0   
## [112] vegan_2.5-7          reprex_1.0.0         digest_0.6.27       
## [115] stats4_4.0.4         munsell_0.5.0        beeswarm_0.2.3      
## [118] viridisLite_0.3.0    vipor_0.4.5          bslib_0.2.4